function suspension
  close all
  L0=1; k=1.e05; m=2000; g=9.81;
  tspan = [0,5];
  initial_w = [L0,0];
  
  [t_vals,sol_vals] = ode45(@(t,w) eom(t,w,k,L0,m,g),tspan,initial_w);
  plot(t_vals,sol_vals(:,1),'Linewidth',2,'Color','r');
  ylim([0 1]);
  set(gca,'FontSize',14);
  xlabel('Time (s)');
  ylabel('Spring length (m)');
  grid on
  pause(3)
  animatespringmass(t_vals,sol_vals,L0);
end
function dwdt = eom(t,w,k,L0,m,g)
   x = w(1); v = w(2);
   dxdt = v;
   dvdt = -k*(x-L0)/m - g;
   dwdt = [dxdt;dvdt];
end
function animatespringmass(times,sols,L)
figure
 for i = 1:length(times)
    clf
    axis([-L L 0 2*L]);
    hold on
    y1 = sols(i,1);
    plot(0,y1,'ks','MarkerSize',60,'MarkerFaceColor','r');
    x1 = [0,0];
    x2 = [0,y1];
    plotspring(x1,x2,0.125);
    pause(0.05)
 end
end
function plotspring(x1,x2,size)
    length = sqrt(dot(x2-x1,x2-x1));
    xx0 = length*[0;2;2.5;3.5;4.5;5.5;6;8]/8.;
    yy0 = size*[0;0;1;-1;1;-1;0;0];
    z0 = length*[1;0];
    z1 = x2-x1;
    theta = acos(dot(z1,z0)/length^2);
    c = z0(1)*z1(2)-z0(2)*z1(1);
    if (c<0) theta = -theta; end
    xx  = cos(theta)*xx0 - sin(theta)*yy0;
    yy =  sin(theta)*xx0 + cos(theta)*yy0;
    xx = xx + x1(1);
    yy = yy + x1(2);
    plot(xx,yy,'k','LineWidth',2);
end